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Abstract 



We examine a square-lattice nearest-neighbor Ising quantum ferromagnet 
coupled to d-dimensional phonon baths. Using the density-matrix equation, 
we calculate the transition rates between configurations, which determines 
the specific dynamic. Applying the calculated stochastic dynamic in Monte 
Carlo simulations, we measure the lifetimes of the metastable state. As the 
magnetic field approaches \H\/J = 2 at low temperatures, the lifetime pref- 
actor diverges because the transition rates between certain configurations ap- 
proaches zero under these conditions. Near \H\/J = 2 and zero temperature, 
the divergent prefactor shows scaling behavior as a function of the field, tem- 
perature, and the dimension of the phonon baths. With proper scaling, the 
simulation data at different temperatures and for different dimensions of the 
baths collapse well onto two master curves, one for \H\/J > 2 and one for 
\H\/J < 2. 

PACS numbers:64.60.Qb,75.60.Jk,05.50.+q,02.50.Ga 



2 



Typeset using REVT^X 



I. INTRODUCTION 



Compared to the exponentials that occur in the expressions for particular physical quan- 
tities, the associated prefactors are often assumed to be too unremarkable and uninteresting 
to examine, so that the study of prefactors has been ignored in many cases. However, this 
is not always true. In this paper, we explore interesting behavior shown by the prefactors of 
the lifetime of the metastable state of a two-dimensional nearest-neighbor Ising ferromagnet 
interacting with phonon baths. Metastability occurs often in many different systems, rang- 
ing from supercooled fluids and vapors [1,2] to spin systems [3] and quantum liquids [4]. In 
those systems some type of weak noise (for example, thermal noise) can drive the system 
from the metastable state to the stable state across a saddle point at an extremely small 
rate. Therefore, the average waiting time before escape from a metastable state is usually 
extremely long. 

Recently, we examined a square- lattice Ising quantum ferromagnet with a phonon (i.e., 
bosonic) bath attached to each spin, in the presence of a longitudinal magnetic field [5,6]. 
The time evolution of the spin system is determined by the linear coupling of the system with 
the phonon baths, and the dynamic (defined by the transition rates between different spin 
configurations) is calculated using the quantum-mechanical density-matrix equation. The 
resulting dynamic is different from the Glauber dynamic [7] , which can be derived in a similar 
way from coupling with fermionic baths [8]. In the present paper, the calculated transition 
rates are applied in dynamic Monte Carlo simulations in order to measure the lifetimes of 
the metastable state. To measure the lifetimes, the initial configuration is set to be all spins 
up, and a magnetic field which favors all spins down is applied. Then the number of spin-flip 
attempts is measured until the system magnetization reaches zero. At a particular value of 
the magnetic field, \H\/J = 2 where J is the nearest-neighbor exchange coupling constant, 
the present dynamic does not allow transitions between certain configurations, which causes 
the lifetime to diverge as the temperature approaches zero and \H\/J -> 2 ± . Since the energy 
barrier between the stable and metastable states is finite, this means that the prefactor of 
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the lifetime must diverge as \H\/J — > 2 1 * 1 . If the Glauber dynamic is applied to the system 
instead, then the lifetime prefactor has a finite, field-independent value for \H\/ J > 2, which 
is different from that for \H\/J < 2 [9,10]. Consequently, the divergence of the prefactor is 
due to the specific dynamic imposed on the spin system. From the dynamic Monte Carlo 
simulations, we demonstrate that this divergent prefactor reveals a scaling behavior which 
is determined by T d and the ratio (\H\ — 2J)/(fcsT), where T is the absolute temperature, 
ks is Boltzmann's constant, and d is the dimension of the phonon baths. 

Quite recently, Maier and Stein investigated the effects of weak noise on the magneti- 
zation reversal rate in an overdamped one-dimensional classical Ginzburg-Landau model of 
finite length [11]. There is no applied field in their model, so thermal fluctuations drive the 
system from one stable configuration (magnetization near +1) to another (magnetization 
near —1) through a saddle point. They analytically calculated the prefactor of the magneti- 
zation reversal rate, and found that the rate prefactor diverges as the system size L reaches 
7r or 2rr, depending on the boundary conditions. A similar prefactor divergence occurs in a 
two-dimensional nonequilibrium model [12]. Shneidman and Nita [13] studied a metastable 
lattice-gas model with nearest-neighbor interactions and continuous-time Metropolis dynam- 
ics. They calculated analytically the prefactor of the lifetime of the metastable state beyond 
the zero-temperature limit. The lifetime prefactor exhibited distinctive peaks as a function 
of field, which disappear at zero temperature. 

The remainder of this paper is organized as follows. In Sec. II we describe our dynamic 
quantum model, and in Sec. Ill we show how to calculate the prefactor analytically and 
numerically. In Sec. IV we present our scaling analysis of the prefactor, and in Sec. V we 
present a discussion and our conclusions. 

II. DYNAMIC QUANTUM MODEL 

To derive the classical dynamic (i.e., transition rates) from the quantum system, we use 
the Hamiltonian: 



4 



Ti = H sp + H P h + ?4p-ph , (1) 
74p = -JE^- ff E°f. ( 2 ) 

^ph = J2 nuJ Q C \°q > ( 3 ) 



"Wsp-ph — a 



" -(^^)( C t-c 9 )e^ , (4) 



where 7i sp is the spin Hamiltonian, 7i p h is the phonon Hamiltonian, and 7i sp _ p h [14] is the 
Hamiltonian describing the simple linear spin-phonon coupling. The first sum in Eq. (2) 
runs over nearest-neighbor sites only on a two-dimensional square lattice, J(> 0) is the 
ferromagnetic exchange coupling constant, cr? are the z components of Pauli spin operators 
attached to lattice site j (in our notation, their eigenvalues are ±1), and H is a longitudinal 
magnetic field. The index q is the wave vector of a phonon mode, uj q is the angular frequency 
of the phonon mode with wave vector q, and and c g are the corresponding creation and 
annihilation operators. The constant A (its dimension is energy) is the coupling strength 
between the spin system and the phonon heat bath, N is the number of unit cells in the 
system, M is the mass per unit cell, and Rj is the position of site j. Details of this model 
were presented in Refs. [5,6], and thus here we only briefly sketch the main ideas. 

With the given spin Hamiltonian, the dynamic is determined by the generalized master 
equation [15,16]: 

+ 

Im'm = ~ , W m = }^ W km , (5) 

where pit) is the time dependent density matrix of the spin system, to', n, k, and to are 
eigenstates of H sp , p(t) m > m ={m'\p(t)\m) , and Wk m is the transition rate from the mth to 
the fcth eigenstate. Assuming that the correlation time of the heat bath is much smaller 
than the times of interest, we integrate over all degrees of freedom of the bath. Then the 
transition rate from the mth to the fcth eigenstate of Ti sp becomes [5,6] 
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W km = -r^2\(n q + l,k |Hsp-ph| n q ,m)\ 2 (n q \p ph \n q ) 5{E m - E k - hu q ) , (6) 
n q,n g 

where E m and E k are the energy eigenvalues of H sp and E m > E k . Here n q is the average 
occupation number of the phonon mode with q, and p p h is the density matrix of the phonon 
bath. We can calculate the transition rate when E m < E k similarly. Eventually we obtain 



W km = 



er]h d+1 c d + 2 



(E k — E m ) 



(7) 



where d is the dimension of the heat bath, 6 = 2 (2rc) for d — 1, 2 (3), rj = M/a d where a 
is the lattice constant, c is the sound velocity, and (3 = l/k B T. The two major differences 
from the Glauber dynamic are the energy difference term in the numerator and the negative 
sign in the denominator. In the limit T — > 0, the transition rates vanish when E m = E k 
(this can occur for \H\ = 2 J, 4J). Despite their unusual form, these transition rates satisfy 
detailed balance. The transition rates for different d scaled by (fcsT) d are shown in Fig. 1(a) 
as functions of (3(E k — E m ). This can be compared with the transition rate for the Glauber 
dynamic, shown in Fig. 1(b). 



III. PREFACTOR: THEORY AND SIMULATION 

We apply the calculated transition rates for the (i-dimensional phonon dynamic in Monte 
Carlo simulations and measure the lifetimes of the metastable state. [Here we consider H < 
only in Eq. (2).] First, we summarize theoretical predictions for the lifetime. As T — > 0, the 
exact prediction [17,18] for the energy barrier between the metastable and stable states is 
given by 

T(H,J) =8J£ c -2\H\(f c -e c +l) , (8) 

where the linear critical droplet size in units of the lattice constant is £ c = \2J/\H\], and 
\x] denotes the smallest integer not less than x. For example, for 2J < \H\ < 4J, £ c = 1, 
and for J < \H\ < 2J, £ c = 2. (At \H\ = 2, the critical droplet size changes.) For £ c > 1, the 
critical droplet is a cluster of overturned spins, which is an £ c x (£ c — 1) rectangle with one 
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additional overturned spin on one of its long sides. For £ c — 1, a single overturned spin is 
the critical droplet. This formula is valid only when 2J/\H\ is not an integer and \H\ < 4J. 
At low temperatures the mean lifetime is the inverse of the probability of escaping from the 
metastable well, which can be written as [17] 

(T)=A(H,T)exp(r/T), (9) 

where A(H,T) is the prefactor. In this equation and hereafter we set ks = 1. We measure 
(r) in units of Monte Carlo spin-flip attempts (MCS). In our case, a unit Monte Carlo step 
is defined to be the inverse of the coefficient of the transition rate [Eq. (7)] multiplied by J d : 
<dr]h d+1 c d+2 /(\ 2 J d ). (For the single-droplet decay mode considered here, (r) is independent 
of the system size if measured in MCS [3,19].) Hereafter, the field H and temperature T are 
given in units of J, and we set J = 1. 

The lifetime can be calculated analytically in terms of the shrinking and growing prob- 
abilities of a droplet, using a rejection-free technique [20] and estimating the first-passage 
time until the system contains a critical droplet and overcome the barrier [9,21]. Then from 
Eqs. (7) and (9), the prefactor A can be calculated [6] as a function of the field with the 
(i-dimensional phonon dynamic, using absorbing Markov chains (AMC) [9,21] in the limit 
T -> 0. For 2 < \H\ < 4, 
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-A) d + (8-5 


\\H\) d 


4(2|i2 


"1 - A) d (% - 2\H\) d 



A(H,T^0,d) = ^ o ,^,J . (10) 

For 1< \H\ < 2, 



\H\' 


* + 2(2- | 


H\) d 


2 d + 3 \H\ d (2 - 


\H\) d 



Near \H\ = 2, both below and above, we find from Eqs. (10) and (11) that for \H\ fixed, 

lim A(H,T,d) = — — y. -t, (12) 

T^O, \H&2 K ' ' 7 d \ \H\ - 2\ d V ' 

where d = 2 d+2 for \H\ > 2 and C x = 2 d+3 for \H\ < 2. On the other hand, as \H\ -> 2 ± 
at nonzero T, with the d- dimensional phonon dynamic, the transition rate, W^m [Eq. (7)], 
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approaches T\ \H\ — 2\ d 1 times a constant. Since the transition rate is inversely proportional 
to the prefactor, we find that for T fixed, 

lim A(H,T,d) = ——^——r , (13) 
|h|-»2±,t>o v ; T\\H\- 2\ d ~ l ' v ; 

where the constant C 2 may depend on d. The exact value of C 2 could not be obtained 
analytically, so we examine it numerically (at the end of Sec. IV). As \H\ — > 2 ± 1 for d = 1 
the prefactor approaches C 2 /T so that there is no field dependence, but for d = 2 and 3 the 
prefactor diverges as C 2 /(T\\H\ - 2|) and C 2 /(T\\H\ - 2| 2 ), respectively. 

To investigate the region in which Eqs. (12) and (13) are applicable, and to examine 
whether there is any scaling behavior in the prefactor near \H\ — 2 and T = 0, we performed 
Monte Carlo simulations with absorbing Markov chains (MCAMC) [9,21] near \H\ = 2 at 
low temperatures, using the <i-dimensional phonon dynamic [Eq. (7)]. Average lifetimes 
were measured over 2000 escapes with the linear system size L = 24 and periodic boundary 
conditions. Since this system size is much larger than the critical droplet size near \H\ = 2, 
the dynamic behavior of the system does not depend on the examined system size. (This 
was confirmed by comparing with simulations for L = 48.) The range of temperatures 
examined was between T = 0.006 and T = 0.10, and the range of fields used was from 
\H\ — 2 = ±10~ 5 to | if | — 2 = ±0.3. At significantly low temperatures in a particular field 
or in a field quite close to | if | = 2 at a particular temperature, the MPFUN package [22] 
was used for arbitrarily high-precision calculations. Equation (9) was used to obtain the 
prefactor from the measured average lifetimes. 

Figure 2 illustrates the analytic and numerical (Monte Carlo simulations at T = 0.1) 
prefactors as functions of the field with the Glauber and the ci-dimensional phonon dynamics. 
The analytic results for the prefactor with the phonon dynamics diverge as |if | — > 2 ± . This 
divergence does not occur with the Glauber dynamic (straight solid lines in Fig. 2). The 
simulation data at T = 0.1 agree well with the analytic results far away from |if| = 2, 
but they start to deviate from the analytic results as |if | approaches 2. This deviation is 
more significant for lower- dimensional phonon baths (smaller d) at a fixed temperature, as 
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shown in Fig. 2. For d = 3 the simulation data at T = 0.1 deviate from the analytic results 
only when the field is much closer to 2 than for the d = 1 simulation data. We also find 
that lower-dimensional phonon dynamics need much lower temperatures to agree with the 
analytic results as \H\ approaches 2. From now on, we concentrate on the phonon dynamics 
only. 

Figure 3 shows log-log plots of the prefactor as a function of field at different temperatures 
and for different d. For d — 1, 2, and 3, the simulation data show two regimes where Eqs. (12) 
and (13) hold, respectively. Roughly speaking (more precise statements will be made in 
Sec. V.), when the field is not too close to \H\ = 2 [about the right halves of Figs. 3 (a), (b), 
and (c)], below (or above) \H\ = 2, the prefactors at different temperatures converge to a 
linear curve. For d — 1, these converging linear curves can be found for larger \\H\ — 2| 
at higher T. The slopes of the converging linear curves (this is easily seen for d = 2 and 
3) are approximately —d, as expected from Eq. (12). When the field is close enough to 
| if | = 2 (the left halves of Fig. 3), the prefactor at fixed temperature again behaves linearly 
with \ \H\ — 2| on a logarithmic scale. As expected from Eq. (13), the slopes of those linear 
curves are approximately — (d — 1) and do not depend on temperature, but the intercepts 
are different for different temperatures. 



IV. SCALING 

Assuming that the prefactor A is a generalized homogeneous function [23] of the field 
and temperature only, we can write it as 

A(X ah h, \ aT T) = X aA A{h, T) , (14) 

where h = \H\ — 2. Choosing A = y^A^, this gives 

Mh,T)=( f ) *±(^), (15) 

where $ ± (a;) are scaling functions for x = h/T ah ^ aT > and x < 0, respectively. Since the 
prefactor has the asymptotic behaviors shown in Eqs. (12) and (13), the scaling functions 



$ ± (a;) should have these properties: 

lim ^(x) = — ^— -. , lim ^(x) = . (16) 

|x|-oo d \x\ d |x|-0 \x\ d ~ l 

From the above properties, the scaling exponents are determined to be oa/ot = d and 
a h /a T = 1, so A(h,T) = (1/T) d <S>*{h/T). 

If we rewrite the prefactor using scaling functions, \l/ ± (x) = |x| d_1( l )± (a:), which have 
these properties 

lim = -i- , lim = C 2 , (17) 

|x|->oo Gi |X| |x|->0 

then the prefactor can be written as 

A ^^ik*i^- <i8) 

Log-log plots of T|/i| d_ M vs \h\/T at different temperatures and for different d are shown 
in Fig. 4. The simulation data for different temperatures collapse well onto each other for 
each value of d and for \H\ > 2 (or \H\ < 2). Because C\ for a particular d and \H\ < 2 is 
the same as C 1 for d+ 1 and > 2 [Eqs. (12) and (17)], the simulation data in the regime 
of large \h\/T for d and \H\ < 2 fall on the data for d + 1 and > 2. Using the fact 
that for large \h\/T and \H\ > 2 (\H\ < 2), Ci multiplied by 2 for d is the same as C\ for 
d-l- 1, we make log-log plots of T\h\ d ~ x A multiplied by 2 d vs \h\/T for different temperatures 
and for different d (see Fig. 5). Then all the simulation data for different temperatures and 
different d collapse well onto two master curves, one for \H\ > 2, and the other for \H\ < 2. 
The saturation value of 2 d T\h\ d - 1 A as \h\/T -> (2 d C 2 ) is approximately 2/3 (we do not 
know if this value is exactly 2/3 or not), so that C 2 ~ 2/3 • l/2 d and C 2 oc \jC\. Using 
this approximate expression for C 2 and Eqs. (17) and (18), we can confirm that the slopes 
of the log-log plots of 2 d T\h\ d ~ 1 A vs x = \h\/T (Fig. 5) for large x are —1 for both \H\ > 2 
and | if | < 2, and that their intercepts at x = 1 for \H\ > 2 (\H\ < 2) are approximately 
log 10 (1/4) (log 10 (1/8)). 
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V. DISCUSSION AND CONCLUSIONS 



Our results indicate that the divergence of the lifetime prefactor is caused by the phonon 
dynamic we imposed, not by any unusual properties of the saddle points themselves, as is 
the case for Ref. [11]. Referring to Fig. 6, we discuss saddle points and the most probable 
paths to escape from the metastable state for different fields. 

For \H\ > 2, the saddle point corresponds to the configuration SA, only one overturned 
spin, so the most probable path to escape from the metastable well is A— >-SA. As \H\ — > 2 + , 
the energy difference between SA and C approaches zero (Fig. 6), so the transition rate 
between them becomes zero. Consequently, the attempt frequency to grow the droplet 
along the downhill path SA^C decreases. This is the reason that the lifetime prefactor 
diverges as \H\ — > 2 + . Then, a new probable path towards SC through SA begins to be 
preferred to the downhill path. At \H\ = 2, the path SA^C is forbidden, and the most 
probable path is directed towards the saddle point SC. 

For | if | < 2, the saddle point corresponds to the configuration SB, an L-shaped cluster of 
three overturned spins. For this field, the most probable path to escape from the metastable 
well is A-^SA^C-^SB. As \H\ — > 2~, the energy difference between SA and C and the 
energy difference between C and SB both approach zero, so that the transition rates between 
them become zero. As a result, the attempt frequency along the most probable path towards 
the saddle point SB decreases, which leads to the prefactor divergence as \H\ — > 2~. Then, 
a new probable path towards SC through SA develops (Fig. 6). Thus, SC starts to compete 
with SB. At | if | = 2, the attempt frequency towards SB is zero, and the most probable path 
is towards the saddle point SC. 

The previously shown scaling behavior of the prefactor can be observed from the simu- 
lation data in Fig. 3. In the log-log plots of A vs \\H\ — 2|, a slope of — (d — 1) is found for 
the regime when — 2|/T — » 0, while a slope of — d is found when — 2\/T — > oo. The 
variable \\H\ — 2|/T (not simply either ||if| — 2| or T) determines how those two regimes 
change with field and temperature. The slope of — (d — 1) occurs at smaller \\H\ — 2| at lower 
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temperatures, in order to reduce \\H\ — 2|/T, compared to higher-temperature data. The 
slope of — d occurs at larger | \H\ — 2| at higher temperatures, in order to increase \ \H\—2\/T, 
compared to lower-temperature data. 

In summary, we have examined the prefactor of the lifetime of the metastable state for 
a square-lattice Ising ferromagnet with a d-dimensional phonon bath attached to each spin. 
Using the transition rates calculated from the density-matrix equation, we demonstrated 
analytically and numerically that the metastable lifetime prefactor diverges as \H\ — > 2, and 
that it also scales as a function of \ \H\ — 2|/T for each value of d, near \H\ = 2 and T = 0. 
This scaling is a result of the fact that the prefactor is a generalized homogeneous function 
[23] of the field and temperature. The divergence and scaling of the prefactor are due to the 
chosen dynamic, not to any non-smooth energy landscape at the particular magnetic field. 
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FIGURES 
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FIG. 1. (a) Transition rate W divided by T d vs the energy difference AE scaled by T for the 
(i-dimensional phonon dynamics. (We use units such as ks = 1-) Here we ignore the proportionality 
constant in W, setting X 2 / @rjh d+1 c d+2 = 1. (b) Transition rate W vs AE/T for the Glauber 
dynamic. 
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FIG. 2. The prefactor A vs field \H\ for the Glauber dynamic and the d-dimensional phonon 
dynamics. The open symbols are from the Monte Carlo simulations at T = 0.1 with the phonon 
dynamics. The curves are from the low-temperature analytic results Eqs. (10) and (11) with 
the phonon dynamics, while the straight solid lines are for the Glauber dynamic. For the d = 1 
simulation data with the phonon dynamic, a lower temperature or a field farther away from \H\ = 2 
is needed to agree with the analytic results, compared to the simulation data for higher d. At 
\H\ = 2, the critical droplet size changes, so the theoretical prediction Eq. (8) is not valid. The 
filled symbols are from the Monte Carlo simulations with the phonon dynamics at \H\ = 2, except 
for the filled triangle, which is from the Glauber dynamic. Here the prefactors at \H\ =2 for both 
the phonon and Glauber dynamics were obtained by extrapolating low-temperature simulation 
data (from T/J = 0.04 to T/J = 0.2) to zero temperature. 
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FIG. 3. Log-log plots of the prefactor A(H,T,d) vs \\H\ — 2| at different temperatures for (a) 
d = 1, (b) d = 2, and (c) d = 3 dimensional phonon dynamics, as the field \H\ approaches 2 
from above (empty symbols) and below (filled symbols). In (c), only data for two temperatures 
are shown for clarity, while the simulations were performed for T = 0.01,0.02,0.03,0.06, and 0.1. 
There are two regimes in these plots: (1) As \\H\ — 2\/T — > (the left-hand side of the plots), the 
slopes of the linear curves for different temperatures become — (d — 1). (2) As \\H\ — 2\/T — > oo 
(the right-hand side of the plots) , the curves for different temperatures collapse onto single linear 
curves whose slopes are —d. These results are predicted from Eqs. (12) and (13). 
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FIG. 4. Log-log plots of T\\H\ - 2\ d ~ 1 A(H, T, d) vs x = \\H\ - 2\/T for different tempera- 
tures (for d = 1, T = 0.006,0.01,0.02,0.03; for d = 2, T = 0.02,0.03,0.04,0.06; and for d = 3, 
T = 0.01, 0.02, 0.03, 0.06, 0.1) and different-dimensional phonon dynamics. For each value of d, the 
upper curves are for \H\ > 2, while the lower curves are for \H\ < 2. The simulation data for 
different temperatures collapse well for each value of d. Also, for large x, the curves for d and 
\H\ < 2 coincide with those for d + 1 and \H\ > 2, and the curves multiplied by 2 for d + 1 and 
> 2 (|iJ| < 2) coincide with the curves for d and > 2 < 2). 
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FIG. 5. Log-log plots of T\\H\ - 2\ d ~ l A(H, T, d) multiplied by 2 d vs \\H\ - 2|/T for different T 
and d. For \H\ > 2 (or \H\ < 2), all the simulation data for different T and d collapse onto a single 
curve. The scaling function for \H\ > 2 is different from that for \H\ < 2. Here the saturation 
value of 2 d T\\H\ - 2\ d ~ 1 A(H,T) as \\H\ - 2[/T -> (that is, 2 d C 2 ) is approximately 2/3, so that 
C 2 w 2/3 • l/2 d . 
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FIG. 6. Schematic diagram that shows the most probable paths of the system for different fields. 
Each box represents one overturned spin. Each configuration is labeled. The configuration labeled 
SA is the saddle point for \H\ > 2, SB for \H\ < 2, and SC for \H\ = 2. The energy difference 
between two configurations is given on the right-hand side. As shown, the energy difference between 
SA and C is zero at \H\ = 2, and so is the energy difference between C and SB. Therefore, at 
\H\ = 2, the path from SA to SB through C is forbidden according to the phonon dynamic, 
Eq. (7). Consequently, the system chooses the saddle point SC to reach the configuration SB 
(A^SA^SC^SB). Then the energy barrier at \H\ = 2 is E sc - E A = 2(8 - 2\H\) = 8 since the 
path from SC to SB is downhill. 
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